!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! the program is meant to design the wall of a supersonic wind 
! tunnel, which provides a mean flow at exit with Mach number
! of Md.
! ref to chapter 16 of 'Gas Dynamics' by Maurice J. Zucrow & 
! Joe D. Hoffmann
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! coded by : B. G.
! created  : 2015-07-21
! revised  : 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! called subroutine:
!   VariableDef (module)
!   Thermo
!   Transonic
!   InterP
!   SymmetryP
!   InverseW
!   Exitline
!    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!   k~j      interger step variables
!   thetai   theta step used in InverseW section
!   interc   interpolation coefficient used to calculate Md c-line
!   atemp    used in Thermo subroutine
!   x~M      these array stored the point paramenters
!   alpd     Mach angle of design Mach number
!   stepd    step length along design Mach c-line
!   mdotk1   mdot through between k1-th Mach line point and x-axis
!   dmdot    mass flow rate between k2-th point and x-axis
!   arou-ay  average paramenters between the point and the front point on the same c-line
!   cmdot    mdot below front point
!   mdot1    mass flow rate through the last specific region c-line
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning :
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program SUPERT_IRRO()
   use VariableDef
   implicit none
!   temperory variables
   integer,allocatable::ffn(:)
   integer::k1,k2,k3,k4,k,i,j
   real*8::thetai,interc,atemp,Mtemp,alptemp,ytemp,alpd,stepd
   real*8::arou,aV,atheta,aalpha,ay,cmdot,mdot1,mdotk1,dmdot
   real*8,allocatable::xbarma(:)
!   Define Output Variable
   real*8,allocatable::x(:,:),y(:,:),Vu(:,:),Vv(:,:),V(:,:),p(:,:),T(:,:),rou(:,:),M(:,:)
   
!********************************** program start *****************************************!
!****************    Initiate variable     ***
! open an file and write header
   open(unit=110,file='F:\Nozzle\Modeling\SUPERT_IRRO\temp.txt')
   write(110,*) 'SUPERT_IRRO Data'
   write(110,*)
   close(110)
! initiate the variables
   el=0.0001
   ev=0.1
   pi=3.141592653
   delta=1
 
   Tt=3000
   Pt=7000000
   g=1.2
   Rg=320.0
   
   yt=1.0 !0.1 !0.025 !1.0
   rtu=2.0 !0.2 !0.050 !2.0
   rtd=0.5 !0.05
   Md=3.0
   thetaa=20*pi/180
   stepd=yt*0.2
   
   wallflag=1
   flagy40=0
   mdotk1=0.0
   
   Ny=11
   Nt=21
   Ni=99
   Nj=99
! allocate memory for array
   allocate(xbarma(Ny))
   allocate(ffn(Ni))
   
   allocate(x(Ni,Nj))
   allocate(y(Ni,Nj))
   allocate(Vu(Ni,Nj))
   allocate(Vv(Ni,Nj))
   allocate(V(Ni,Nj))
   allocate(p(Ni,Nj))
   allocate(T(Ni,Nj))
   allocate(rou(Ni,Nj))
   allocate(M(Ni,Nj))
   
!***********************  initial line ******
   call Transonic(xbarma,x(1:Ny,1),y(1:Ny,1),Vu(1:Ny,1),Vv(1:Ny,1),V(1:Ny,1)&
         &,M(1:Ny,1),T(1:Ny,1),p(1:Ny,1),rou(1:Ny,1))
   ffn(1)=1

!**********************   initial region *********
   do k1=2,11
      do k2=2,2*k1-2
         x1=x(k1,k2-1)
         y1=y(k1,k2-1)
         Vu1=Vu(k1,k2-1)
         Vv1=Vv(k1,k2-1)
         x2=x(k1-1,k2-1)
         y2=y(k1-1,k2-1)
         Vu2=Vu(k1-1,k2-1)
         Vv2=Vv(k1-1,k2-1)
         
         call InterP()
         if (flagy40==1) then
            flagy40=0
            exit
         end if
         
         x(k1,k2)=x4
         y(k1,k2)=y4
         Vu(k1,k2)=Vu4
         Vv(k1,k2)=Vv4
         V(k1,k2)=V4
         M(k1,k2)=M4
         T(k1,k2)=T4
         p(k1,k2)=p4
         rou(k1,k2)=rou4
      end do
      x1=x(k1,k2-1)
      y1=y(k1,k2-1)
      Vu1=Vu(k1,k2-1)
      Vv1=Vv(k1,k2-1)
      
      call SymmetryP()
      
      x(k1,k2)=x4
      y(k1,k2)=y4
      Vu(k1,k2)=Vu4
      Vv(k1,k2)=Vv4
      V(k1,k2)=V4
      M(k1,k2)=M4
      T(k1,k2)=T4
      p(k1,k2)=p4
      rou(k1,k2)=rou4
      ffn(k1)=k2
   end do
!   check mass flow rate balance
   mdot1=0.0
   do i=2,2*Ny-1
      arou=(rou(k1-1,i-1)+rou(k1-1,i))/2
      aV=(V(k1-1,i-1)+V(k1-1,i))/2
      atheta=(atan(Vv(k1-1,i-1)/Vu(k1-1,i-1))+atan(Vv(k1-1,i)/Vu(k1-1,i)))/2
      aalpha=(asin(1/M(k1-1,i-1))+asin(1/M(k1-1,i)))/2
      ay=(y(k1-1,i-1)+y(k1-1,i))/2
      mdot1=mdot1+2*pi*ay*arou*aV*sin(aalpha)/sin(aalpha-atheta)*(y(k1-1,i-1)-y(k1-1,i))
   end do
   
! write the initial region data
   open(unit=110,file='F:\Nozzle\Modeling\SUPERT_IRRO\temp.txt',access='append',status='old')
   write(110,*) 'Mass flow rate for initial line, mdot mdot1d CDm'
   write(110,*) mdot,mdot1d,CDm
   write(110,*)
   write(110,*) 'Initial Region Points, x y Vu Vv V M T p rou'
   write(110,*)
   do i=1,Ny
      do j=1,2*i-1
         write(110,*) i,j,x(i,j),y(i,j),Vu(i,j),Vv(i,j),V(i,j),M(i,j),T(i,j),p(i,j),rou(i,j)
      end do
   end do
   write(110,*)
   write(110,*) 'Inverse Wall Region Points, i j x y Vu Vv V M T p rou'
   write(110,*)
   close(110)

!****************************** inverse wall region   ************
   thetai=thetaa/(Nt-1)   ! initiate the step angle
   theta4=0.0   ! initiate the wall angle
   k=2*Ny-1   ! initiate the first inverse wall c-line point number
   do k1=Ny+1,Ny+Nt-1
      theta4=theta4+thetai	!get the wall angle used in the cycle
      ! there the next c+ characteristics line must not intersect the wall before the prescribed point
      do k3=2,k
         x(k1,1)=rtd*sin(theta4)
         y(k1,1)=yt+rtd*(1-cos(theta4))
         
         x3=x(k1-1,1)
         y3=y(k1-1,1)
         Vu3=Vu(k1-1,1)
         Vv3=Vv(k1-1,1)
         x4=x(k1,1)
         y4=y(k1,1)
         
         call Thermo(V(k1-1,k3),p(k1,1),rou(k1,1),T(k1,1),atemp,Mtemp)
         alptemp=atan(Vv(k1-1,k3)/Vu(k1-1,k3))+asin(1/Mtemp)
         ytemp=y(k1-1,k3)+(x(k1,1)-x(k1-1,k3))*tan(alptemp)
         
         if (ytemp < y(k1,1)) then
            x1=x(k1-1,k3)
            y1=y(k1-1,k3)
            Vu1=Vu(k1-1,k3)
            Vv1=Vv(k1-1,k3)
            
            call InverseW()
            
            Vu(k1,1)=Vu4
            Vv(k1,1)=Vv4
            V(k1,1)=V4
            M(k1,1)=M4
            T(k1,1)=T4
            p(k1,1)=p4
            rou(k1,1)=rou4
            exit
         end if
      end do
      
!     successive points
      do k2=2,k-k3+2
         x1=x(k1,k2-1)
         y1=y(k1,k2-1)
         Vu1=Vu(k1,k2-1)
         Vv1=Vv(k1,k2-1)
         x2=x(k1-1,k2+k3-2)
         y2=y(k1-1,k2+k3-2)
         Vu2=Vu(k1-1,k2+k3-2)
         Vv2=Vv(k1-1,k2+k3-2)
         
         call InterP()
         if (flagy40==1) then
            flagy40=0
            exit
         end if
         
         x(k1,k2)=x4
         y(k1,k2)=y4
         Vu(k1,k2)=Vu4
         Vv(k1,k2)=Vv4
         V(k1,k2)=V4
         M(k1,k2)=M4
         T(k1,k2)=T4
         p(k1,k2)=p4
         rou(k1,k2)=rou4
      end do
      k=k2
!      axial point
      x1=x(k1,k2-1)
      y1=y(k1,k2-1)
      Vu1=Vu(k1,k2-1)
      Vv1=Vv(k1,k2-1)
      
      call SymmetryP()
      
      x(k1,k2)=x4
      y(k1,k2)=y4
      Vu(k1,k2)=Vu4
      Vv(k1,k2)=Vv4
      V(k1,k2)=V4
      M(k1,k2)=M4
      T(k1,k2)=T4
      p(k1,k2)=p4
      rou(k1,k2)=rou4
      ffn(k1)=k2
      
!      if the Md attained, stop and exit the loop
      if (M4 >= Md) exit

   end do
   k4=k1   ! store the ending c-line number
!   if the Md can not calculated from the region provided circular arc angle
   if (M4 < Md) then
      write(*,*) 'The wall slope at the end point of circular arc too small !'
      stop
   end if
!   calculate Md chracteristics line and store it as the last c-line in the region
   interc=(Md-M(k1-1,ffn(k1-1)))/(M(k1,ffn(k1))-M(k1-1,ffn(k1-1)))   !interpolate coefficient
   theta4=theta4-thetai+thetai*interc
!   first point, use circle function to insure the calculated point in the arc
   x(k1,1)=rtd*sin(theta4)
   y(k1,1)=yt+rtd*(1-cos(theta4))
   Vu(k1,1)=Vu(k1-1,1)+(Vu(k1,1)-Vu(k1-1,1))*interc
   Vv(k1,1)=Vv(k1-1,1)+(Vv(k1,1)-Vv(k1-1,1))*interc
   V(k1,1)=sqrt(Vu(k1,1)**2+Vv(k1,1)**2)
   call Thermo(V(k1,1),p(k1,1),rou(k1,1),T(k1,1),atemp,M(k1,1))
   mdot1=0.0   ! mdot one- check mass flow rate balance
!   successive points
   do k2=2,k-1
      x(k1,k2)=x(k1-1,k2+k3-2)+(x(k1,k2)-x(k1-1,k2+k3-2))*interc
      y(k1,k2)=y(k1-1,k2+k3-2)+(y(k1,k2)-y(k1-1,k2+k3-2))*interc
      Vu(k1,k2)=Vu(k1-1,k2+k3-2)+(Vu(k1,k2)-Vu(k1-1,k2+k3-2))*interc
      Vv(k1,k2)=Vv(k1-1,k2+k3-2)+(Vv(k1,k2)-Vv(k1-1,k2+k3-2))*interc
      V(k1,k2)=sqrt(Vu(k1,k2)**2+Vv(k1,k2)**2)
      call Thermo(V(k1,k2),p(k1,k2),rou(k1,k2),T(k1,k2),atemp,M(k1,k2))
      arou=(rou(k1,k2-1)+rou(k1,k2))/2
      aV=(V(k1,k2-1)+V(k1,k2))/2
      atheta=(atan(Vv(k1,k2-1)/Vu(k1,k2-1))+atan(Vv(k1,k2)/Vu(k1,k2)))/2
      aalpha=(asin(1/M(k1,k2-1))+asin(1/M(k1,k2)))/2
      ay=(y(k1,k2-1)+y(k1,k2))/2
      mdot1=mdot1+2*pi*ay*arou*aV*sin(aalpha)/sin(aalpha-atheta)*(y(k1,k2-1)-y(k1,k2))
   end do
!   last point
   x(k1,k2)=x(k1-1,k2+k3-3)+(x(k1,k2)-x(k1-1,k2+k3-3))*interc
   y(k1,k2)=y(k1-1,k2+k3-3)+(y(k1,k2)-y(k1-1,k2+k3-3))*interc
   Vu(k1,k2)=Vu(k1-1,k2+k3-3)+(Vu(k1,k2)-Vu(k1-1,k2+k3-3))*interc
   Vv(k1,k2)=Vv(k1-1,k2+k3-3)+(Vv(k1,k2)-Vv(k1-1,k2+k3-3))*interc
   V(k1,k2)=sqrt(Vu(k1,k2)**2+Vv(k1,k2)**2)
   call Thermo(V(k1,k2),p(k1,k2),rou(k1,k2),T(k1,k2),atemp,M(k1,k2))
   arou=(rou(k1,k2-1)+rou(k1,k2))/2
   aV=(V(k1,k2-1)+V(k1,k2))/2
   atheta=(atan(Vv(k1,k2-1)/Vu(k1,k2-1))+atan(Vv(k1,k2)/Vu(k1,k2)))/2
   aalpha=(asin(1/M(k1,k2-1))+asin(1/M(k1,k2)))/2
   ay=(y(k1,k2-1)+y(k1,k2))/2
   mdot1=mdot1+2*pi*ay*arou*aV*sin(aalpha)/sin(aalpha-atheta)*(y(k1,k2-1)-y(k1,k2))

   xa=x(k1,1)
   ya=y(k1,1)
   thetaa=atan(Vv(k1,k2)/Vu(k1,k2))
! write the data in the file
   open(unit=110,file='F:\Nozzle\Modeling\SUPERT_IRRO\temp.txt',access='append',status='old')
   do i=Ny+1,k4
      do j=1,ffn(i)
         write(110,*) i,j,x(i,j),y(i,j),Vu(i,j),Vv(i,j),V(i,j),M(i,j),T(i,j),p(i,j),rou(i,j)
      end do
   end do
   write(110,*)
   write(110,*) 'Design Wall Region Points, i j x y Vu Vv V M T p rou'
   write(110,*)
   close(110)
   
!*************************** calculate the region R and the wall location   ***************
   alpd=asin(1/M(k1,k2))
   do k1=k4+1,Ni
!   first point calculation
      x2=x(k1-1,ffn(k1-1))
      y2=y(k1-1,ffn(k1-1))
      Vu2=Vu(k1-1,ffn(k1-1))
      Vv2=Vv(k1-1,ffn(k1-1))
      
      call Exitline(alpd,stepd)
      
      x(k1,ffn(k1-1))=x4
      y(k1,ffn(k1-1))=y4
      Vu(k1,ffn(k1-1))=Vu4
      Vv(k1,ffn(k1-1))=Vv4
      V(k1,ffn(k1-1))=V4
      M(k1,ffn(k1-1))=M4
      T(k1,ffn(k1-1))=T4
      p(k1,ffn(k1-1))=p4
      rou(k1,ffn(k1-1))=rou4
      
      arou=(rou(k1,ffn(k1-1))+rou(k1-1,ffn(k1-1)))/2
      aV=(V(k1,ffn(k1-1))+V(k1-1,ffn(k1-1)))/2
      atheta=(atan(Vv(k1,ffn(k1-1))/Vu(k1,ffn(k1-1)))+atan(Vv(k1-1,ffn(k1-1))/Vu(k1-1,ffn(k1-1))))/2
      aalpha=(asin(1/M(k1,ffn(k1-1)))+asin(1/M(k1-1,ffn(k1-1))))/2
      ay=(y(k1,ffn(k1-1))+y(k1-1,ffn(k1-1)))/2
      cmdot=mdotk1	! record mass flow rate on the exitline before the previous point
!      mass flow rate on the exitline before the current point
      mdot1=mdot1+2*pi*ay*arou*aV*sin(aalpha)/sin(aalpha+atheta)*(y(k1,ffn(k1-1))-y(k1-1,ffn(k1-1)))
      dmdot=mdotk1
!         if it arrived the end point, that means y2<y1
      if (dmdot>=mdot) then
         ffn(k1)=1
         k2=ffn(k1-1)
         interc=(mdot-cmdot)/(dmdot-cmdot)
         x(k1,1)=x(k1-1,k2)+(x(k1,k2)-x(k1-1,k2))*interc
         y(k1,1)=y(k1-1,k2)+(y(k1,k2)-y(k1-1,k2))*interc
         Vu(k1,1)=Vu(k1-1,k2)+(Vu(k1,k2)-Vu(k1-1,k2))*interc
         Vv(k1,1)=Vv(k1-1,k2)+(Vv(k1,k2)-Vv(k1-1,k2))*interc
         V(k1,1)=sqrt(Vu(k1,k2)**2+Vv(k1,k2)**2)
         call Thermo(V(k1,1),p(k1,1),rou(k1,1),T(k1,1),atemp,M(k1,1))
         exit
      end if
!      successive ponts
      do k2=ffn(k1-1)-1,1,-1
         x1=x(k1,k2+1)
         y1=y(k1,k2+1)
         Vu1=Vu(k1,k2+1)
         Vv1=Vv(k1,k2+1)
         x2=x(k1-1,k2)
         y2=y(k1-1,k2)
         Vu2=Vu(k1-1,k2)
         Vv2=Vv(k1-1,k2)
         
         call InterP()
         if (flagy40==1) then
            flagy40=0
            exit
         end if
        
         x(k1,k2)=x4
         y(k1,k2)=y4
         Vu(k1,k2)=Vu4
         Vv(k1,k2)=Vv4
         V(k1,k2)=V4
         M(k1,k2)=M4
         T(k1,k2)=T4
         p(k1,k2)=p4
         rou(k1,k2)=rou4
         
         arou=(rou(k1,k2)+rou(k1,k2+1))/2
         aV=(V(k1,k2)+V(k1,k2+1))/2
         atheta=(atan(Vv(k1,k2)/Vu(k1,k2))+atan(Vv(k1,k2+1)/Vu(k1,k2+1)))/2
         aalpha=(asin(1/M(k1,k2))+asin(1/M(k1,k2+1)))/2
         ay=(y(k1,k2)+y(k1,k2+1))/2
         cmdot=dmdot   ! mass flow rate before previous point
!         mass flow rate before current point
         mdot1=mdot1+2*pi*ay*arou*aV*sin(aalpha)/sin(aalpha-atheta)*(y(k1,k2)-y(k1,k2+1))
         if (dmdot>=mdot) then
            ffn(k1)=ffn(k1-1)-k2+1
            interc=(mdot-cmdot)/(dmdot-cmdot)
            x(k1,k2)=x(k1,k2+1)+(x(k1,k2)-x(k1,k2+1))*interc
            y(k1,k2)=y(k1,k2+1)+(y(k1,k2)-y(k1,k2+1))*interc
            Vu(k1,k2)=Vu(k1,k2+1)+(Vu(k1,k2)-Vu(k1,k2+1))*interc
            Vv(k1,k2)=Vv(k1,k2+1)+(Vv(k1,k2)-Vv(k1,k2+1))*interc
            V(k1,k2)=sqrt(Vu(k1,k2)**2+Vv(k1,k2)**2)
            call Thermo(V(k1,k2),p(k1,k2),rou(k1,k2),T(k1,k2),atemp,M(k1,k2))
            exit
         end if
      end do
!      data translating
      do i=1,ffn(k1)
         x(k1,i)=x(k1,i+k2-1)
         y(k1,i)=y(k1,i+k2-1)
         Vu(k1,i)=Vu(k1,i+k2-1)
         Vv(k1,i)=Vv(k1,i+k2-1)
         V(k1,i)=V(k1,i+k2-1)
         M(k1,i)=M(k1,i+k2-1)
         T(k1,i)=T(k1,i+k2-1)
         p(k1,i)=p(k1,i+k2-1)
         rou(k1,i)=rou(k1,i+k2-1)
      end do
   end do
! write the data in the file
   open(unit=110,file='F:\Nozzle\Modeling\SUPERT_IRRO\temp.txt',access='append',status='old')
   do i=k4+1,k1
      do j=1,ffn(i)
         write(110,*) i,j,x(i,j),y(i,j),Vu(i,j),Vv(i,j),V(i,j),M(i,j),T(i,j),p(i,j),rou(i,j)
      end do
   end do
   close(110)
   write(*,*) 'complete !'
!*************************************program end *********************************************!
end program SUPERT_IRRO
